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Genomic regions displaying outstanding correlation with some environmental variables are likely to be 
under selection and this is the rationale of recent methods of identifying selected loci and retrieve functional 
information about them. To be efficient, such methods need to be able to disentangle the potential effect 
of environmental variables from the confounding effect of population history. For the routine analysis of 
genomewide data-sets, one also need fast inference and model selection algorithms. We describe a method 



ba sed on an explicit spat ial model that builds on the theoretical and computational framework developed 



by 



Rue et al 



(2009) and 



Lindgren et al 



(|201l[ ). The methods allows one to quantify correlation between 



genotypes and environmental variables and to rank loci accordingly. It works for SNP and AFLP data 
obtained either at the individual or at the population level. We provide R scripts with detailed comments 
that can be used readily for the analysis of real data without specific prior knowledge of the R language. 
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Background 



Detecting genomic regions under selection provides important insight about evolutionary processes in humans 
and other organisms. This task has relied exte nsively o n the analysis of either population differentiation, allele 
frequency spectrum or haplotype structure ([Nielsen! 120051) . Alternatively, one can look at loci displaying 
outstanding correlation with some environmental variables. This direct approach has the advantage to 
provide important functional information. However, developing an accurate and fast method toward this goal 
offers a number of statistical challenges: (i) the method should be able to deal with individual or population 
data, any kind of markers (prominently SNP and AFLP data) and any kind of quantitative or categorical 



environmental variables; ( ii) the method should not rely on any ad- hoc method (such as the Mantel test, 



Guillot and Rousset 



20121 ) and be based as much as possible on a well grounded statisti cal model; (hi) the 



method should account for spat ial auto-correlation patterns due to population history (jGuillot et al 



Novembre and Di Rienzo 



2009; 



20091 ); (iv) the method should be fast. 



Joost et al.l (120071 120081 ) proposed a method based on a logistic regression that addresses points (i-ii) 



and ICoop et al.l (|2010l ) proposed an improvement by tackling poin t (iii) with a hi e rarch ical model within an 



MCMC framework. The latter method has been implemented by 



Hancock et al 



honk 



h008i) and 



Hancock et al 



We describe here a spatially explicit and MCMC-free method that addresses points (i-iii) but also solves the 
computing time issue (iv) which, at the era of genome-wide data-sets, remains a crucial bottleneck. 

Method proposed 



We denote by (zi)i = i / a collection of allele counts obtained at / geographical locations (with haploid 

sample sizes nt) and by (j/i)»=i / some measurements of an environmental variable obtained at the same 



set of geographical locations. Denoting by fa the local allele frequency at geographical location Sj, we model 



the dependency of / and y as 

fi = logit~ 1 (aa; l + by, + c). (1) 

where x is an unobserved spatially correlated random effect that accounts for spatial auto-correlation due 
to population history. This random effect x is assumed to be a Gaussian random field, i.e. at any finite 
collection of sites, the entries (x%, ...,x n ) form jointly a Gaussian random vector. The distribution of x is 
characterised by its mean function /i(s) = E\x(s)] and it spatia l covariance function c(s, s') — Cov[x(s), x(s')}. 



As commonly done in spatial statistics (jGelfand et all l2010f) , we make the assumption that x is isotropic 



and stationary, namely that /x(s) is constant and c(s, s') depends only on the spatial lag h = \s — s'|. Further, 
we assume that c(s, s') belongs to the Matern family i.e. 



c(h) 



-{nh) v K v {vh) 



(2) 



where K v is the modified Bessel function of the second kind and order i/>0,K>0isa scaling parameter 
and a 2 is the marginal va riance. For computational reasons the smoothness parameter v is taken equal to 1 



(see 



Lindgren et al 



2011 



for details). 



For co-dominant markers, and assuming that Hardy- Weinberg equilibrium holds locally, the allele counts 
are sampled from a binomial distribution: 



Zi ~ Binom(ni, fi) 



(3) 



For bi- allelic dominant markers such as AFLP markers, 



Zi - Binom(?ij,/-) 



(4) 



where f[ — 2/j(l — fi) + fi - The model proposed appears t herefore as a speci al case of a spatial generalized 



(|2007l ). the latter being obtained by 



linear mixed model (SGLMM). It generalizes the model of Uoost et al 
dropping x in Eq. Q] 

A key property is that under the weak assumption that x is a stationary Gaussian random field with 

a Matern covariance function, inference under the SGLM model above (Eq. [l]|4|) can be performed within 
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the theoretical and computational framework developed by IRue et al.l (|2009D and ILindgren et al.l ([201 1[ ) . 



Briefly, the former develops a framework for Bayesian inference in a broad class of models enjoying a latent 
Gaussian structure. The latter, bridges a gap between Markov random fields and Gaussian random fields 
making it possible to combine flexibility for modelling of Gaussian random fields and efficiency of Markov 
random fields for inference. Following the phrasing of these authors, the combined approach will be referred 
below to as SPDE-INLA. Casting the present problem in the framework of SPDE-1NLA opens the door to 
accurate and extremely fast methods of detecting selection using the r-inla package. See references above 
and www . r-inla . org for details. 

Denoting by Tr(.\m) the probability distribution under some model m and 9 a vector of unknown param- 
eters, the r-inla package provides an evaluation of the integrated likelihood 



j(m) = n(z\m) = J ir(z\9,m)ir(6\m)d6 



(5) 



This quantity can be used to rank loci by de creasing evidence of selec tion or to test Hq : b = against 



H i : b ^ using L(mi) / L(mo) as a test statistic (jKass and Rafterv 



19951) as demonstrated by 



Coop et al 



d2010h . 



Discussion 



One limit common to the approach proposed here and that of lJoost et alj (|20071 ) and lCoop et al.l ([20101 ) is that 



loci are treated separately and assumed to be conditionally independent (no residual linkage disequilibrium 
not accounted for by x and y). This assumption will be clearly violated for dense SNPs data-sets. This 
aspect requires more work for a rigorous and efficient control of false discovery. However we note that 
Bonferonni-type correction offers a solution to protect oneself against false positives. Moreover, potential 
linkage disequilibrium not accounted for has no effect affect on the ranking of the loci in terms of evidence 



of selection. The approach proposed can be therefore readily used to identify conspicuous loci that are likely 



to be the target of selection. 



The approach proposed embeds the main features of the models of ICoop et al.l (|2010j ) and we do not 
expect any major improvement over the latter in terms of accuracy. However, the SPDE-INLA method is 
MCMC-free and extremely fast. For example, computation of Bayesian estimate and integrated likelihood 
can be performed in about 5 seconds per locus on a PC with a dual-core 3GHz chip-set for a data-set 
containing 200 geographical locations. This has to be compared to minutes- to hour-long computing time of 
any MCMC method. The SPDE-INLA method is therefore a promising approach for the routine analysis of 
genomewide data-sets. 
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R code for inference and model choice with SNP data 



The R code belows assumes that the R-inla package is installed. See www.r-inla.org for instructions. 
The main task achieved is to compute the marginal likelihood in four sub-models. The (simulated) dataset 
analyzed here consists of allele counts at a bi-allelic locus for 500 individuals located at 200 geographical 
sites. The marker is assumed to be co-dominant. The toy data used here are availbale from 



http : //www2 . imm.dtu.dk/- gigu/web_inla-spde_selection/. 



############################################################################################# 

library(INLA) 

## loading data 

# geographical sampling location: a 2-column matrix 
s <- as .matrix (read. table ( 'toy-data/s .txt ') ) 

colnames(s) <- c( 'xcoord' , 'ycoord' ) 

# environmental variable: a vector 
y <- as .matrix (read. table ( 'toy-data/y .txt ' ) ) 

colnames(y) <- 'y' 

# allele counts at some SNP markers: a vector 
z <- as. matrix (read.table('toy-data/z-codom.txt')) 

colnames(z) <- 'z' 

# haploid population size at the various sites: a vector 
popsize <- as . vector (read. table( 'toy-data/pop-size .txt ') ) 

colnames (popsize) <- 'popsize' 

## creating mesh suitable for inal-spde computations 
mesh <- inla. mesh. create . helper (points=s , 

off set=c(0. 1, .5), 

max. edge=c (0.06,1)) 

spde = inla. spde . create (mesh, model='matern' , param=list (alpha=2) ) 
data.df <- data. frame (xcoord=s [, 1] , 

ycoord=s [,2] , 

z=z , 

Ntrials=popsize , 

y=y> 

spatial=mesh$idx$loc) 

## 

## Compute parameter estimates and integrated likelihood 

## model z~l 

formula = z ~ 1 

res <- inla(f ormula, 

f amily=' binomial' , 

Ntrials=popsize , 

data = data.df, 

control. data = list (link= ' logit ' ) , 
control . fixed= list (prec . intercept = 1), 
control . compute = list(hyperpar=FALSE, 

return . marginal s=FALSE , ml ik=TRUE) , 
control . inla=list (tolerance=tolerance) , 
verbose=TRUE) 
lmlik.l <- res$mlik[l,l] 
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## model z~l+y 
formula = z ~ 1 + y 
res <- inla(f ormula, 

f amily=' binomial' , 

Ntrials=popsize , 

data = data.df , 

control . data = list (link= ' logit ' ) , 
control . fixed= list (prec . intercept = 1) , 
control . compute = list(hyperpar=FALSE, 

return . marginal s=FALSE , ml ik=TRUE) , 
control . inla=list (tolerance=tolerance) , 
verbose=TRUE) 

lmlik.ly <- res$mlik [1 , 1] 

## estimate of b parameter 

b.est.ly <- res$summary. fixed ['y' , 'mean'] 

## model z~l+x 

formula = z ~ 1 + f (spatial ,model=spde) 
res <- inla(f ormula, 

f amily=' binomial' , 

Ntrials=popsize , 

data = data.df, 

control. data = list (link= ' logit ') , 
control . fixed= list (prec . intercept = 1) , 
control . compute = list(hyperpar=FALSE, 

return . marginal s=FALSE , ml ik=TRUE) , 
control . inla=list (tolerance=tolerance) , 
verbose=TRUE) 
lmlik.lx <- res$mlik[l,l] 

## model z~l+x+y 

formula = z ~ 1 + f (spatial ,model=spde) + y 
res <- inla(f ormula, 

f amily=' binomial' , 

Ntrials=popsize , 

data = data.df, 

control. data = list (link= ' logit ') , 
control. fixed= list (prec . intercept = 1) , 
control . compute = list(hyperpar=FALSE, 

return . marginal s=FALSE , ml ik=TRUE) , 
control . inla=list (tolerance=tolerance) , 
verbose=TRUE) 

lmlik.lxy <- res$mlik [1 , 1] 

## estimate of b parameter 

b.est.lxy <- res$summary .f ixed['y' , 'mean'] 
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R code for inference and model choice with AFLP data 



The (simulated) dataset analyzed here consists of allele counts at a single bi-allelic locus for 200 individuals 
located at 200 geographical sites. The marker is assumed to be dominant. 



############################################################################################# 
## loading data 

# geographical sampling location: a 2-column matrix 
s <- as .matrix (read. table ( 'toy-data/s .txt ') ) 

colnames(s) <- c( 'xcoord' , 'ycoord' ) 

# environmental variable : a vector 
y <- as. matrix (read. table ( 'toy-data/y. txt ') ) 

colnames(y) <- 'y' 

# allele counts at some SNP markers: a vector 
z <- as .matrix (read. table ( 'toy-data/z-dom. txt ') ) 

colnames(z) <- 'z' 

# haploid population size at the various sites: a vector 

ploidy <- 2 



## creating mesh suitable for inal-spde computations 
mesh <- inla. mesh. create. helper (points=s, 

off set=c(0. 1, .5), 

max. edge=c (0.06,1)) 

spde = inla. spde . create (mesh, model="matern" , param=list (alpha=2) ) 
data.df <- data. frame (xcoord=s [, 1] , 

ycoord=s [,2] , 

z=z , 

Ntrials=popsize , 

y=y. 

spatial=mesh$idx$loc) 

## Compute parameter estimates and integrated likelihood 
## under model z~l+x+y 

## note the use of the cbinomial (clumpped binomial) distribution to account for 
## the dominant nature of AFLP markers 
formula = z ~ 1 + f (spatial ,model=spde) + y 
res <- inla(f ormula, 

f amily=" cbinomial" , 

Ntrials=ploidy , 

data = data.df, 

control. data = list(link="logit") , 

control . fixed= list (prec . intercept = 1), 

control . compute = list(hyperpar=FALSE, 
return . marginal s=FALSE , ml ik=TRUE) , 

control . inla=list (tolerance=tolerance) , 

verbose=TRUE) 
lmlik.lxy <- res$mlik [1 , 1] 
## estimate of b parameter 
b.est.lxy <- res$summary .f ixed['y' , 'mean'] 
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